function [ n, result ] = Runge_Kutta( x, h, y0, funct, orde_runge_kutta)
%RUNGE_KUTTA Summary of this function goes here
%   Detailed explanation goes here

n = (x(2) - x(1)) / h;
tempx = x(1);
tempy = y0;

result = 0; % inisialisasi

% Runge Kutta Orde 3
% ------------------------------------------------------------%
if (orde_runge_kutta == 3) 
    % Fungsi 2 variabel
    % ------------------------------------------------------------%
    if (nargin(funct) == 2)
        for r = 1:n
            k1 = h * funct(tempx, tempy);
            k2 = h * funct(tempx + (h/2), tempy + (k1/2));
            k3 = h * funct(tempx + h, tempy - k1 + (2 * k2));
            tempy = tempy + ((k1 + (4 * k2) + k3) / 6);
            tempx = tempx + h;
        end
        result = tempy;
    % Fungsi 3 variabel
    % ------------------------------------------------------------%
    elseif (nargin(funct) == 3)
        tempy = y0(1);
        tempz = y0(2);
        
        for r = 1:n
            k1 = h * funct(tempx, tempy, tempz);
            k2 = h * funct(tempx + (h/2), tempy + (k1/2), tempz + (k1/2));
            k3 = h * funct(tempx + h, tempy - k1 + (2 * k2), tempz - k1 + (2 * k2));
            tempz = tempz + ((k1 + (4 * k2) + k3) / 6);
            tempy = tempy + ((k1 + (4 * k2) + k3) / 6);
            tempx = tempx + h;
        end
        
        result(1) = tempy;
        result(2) = tempz;
    end
% Runge Kutta Orde 4
% ------------------------------------------------------------%
elseif (orde_runge_kutta == 4)
    % Fungsi 2 variabel
    % ------------------------------------------------------------%
    if (nargin(funct) == 2) 
        for r = 1:n
            k1 = h * funct(tempx, tempy);
            k2 = h * funct(tempx + (h/2), tempy + (k1/2));
            k3 = h * funct(tempx + (h/2), tempy + (k2/2));
            k4 = h * funct(tempx + h, tempy + k3);
            tempy = tempy + ((k1 + (2 * k2) + (2 * k3) + k4) / 6);
            tempx = tempx + h;
        end
        result = tempy;
    % Fungsi 3 variabel
    % ------------------------------------------------------------%
    elseif (nargin(funct) == 3) 
        tempy = y0(1);
        tempz = y0(2);
        
        for r = 1:n
            k1 = h * funct(tempx, tempy, tempz);
            k2 = h * funct(tempx + (h/2), tempy + (k1/2), tempz + (k1/2));
            k3 = h * funct(tempx + (h/2), tempy + (k2/2), tempz + (k2/2));
            k4 = h * funct(tempx + h, tempy + k3, tempz + k3);
            tempz = tempz + ((k1 + (2 * k2) + (2 * k3) + k4) / 6);
            tempy = tempy + ((k1 + (2 * k2) + (2 * k3) + k4) / 6);
            tempx = tempx + h;
        end
        
        result(1) = tempy;
        result(2) = tempz;
    end
end 

end
